stanEx3
ex8
non linear regression
with one explanatory variable fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}single linear regressionex8-0.stan
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}n=20
x=runif(n,0,20)
y=rnorm(n,x*2+5,5)
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -8616.18
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 18 -41.7695 0.00198504 0.00147553 0.9751 0.9751 36
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -41.77
## b0 7.04
## b1 1.85
## s 4.90
## m0[1] 29.02
## m0[2] 20.97
## m0[3] 38.71
## m0[4] 15.86
## m0[5] 15.97
## m0[6] 16.76
## m0[7] 37.52
## m0[8] 38.07
## m0[9] 24.05
## m0[10] 21.82
## m0[11] 15.88
## m0[12] 30.32
## m0[13] 27.14
## m0[14] 33.01
## m0[15] 39.34
## m0[16] 8.05
## m0[17] 31.11
## m0[18] 28.46
## m0[19] 30.55
## m0[20] 43.41
## y0[1] 28.69
## y0[2] 15.77
## y0[3] 31.38
## y0[4] 14.15
## y0[5] 26.53
## y0[6] 20.85
## y0[7] 41.17
## y0[8] 42.76
## y0[9] 19.51
## y0[10] 28.79
## y0[11] 15.93
## y0[12] 28.48
## y0[13] 23.36
## y0[14] 38.88
## y0[15] 36.94
## y0[16] 10.18
## y0[17] 33.10
## y0[18] 26.09
## y0[19] 34.13
## y0[20] 33.11
## m1[1] 8.05
## m1[2] 11.98
## m1[3] 15.90
## m1[4] 19.83
## m1[5] 23.76
## m1[6] 27.69
## m1[7] 31.62
## m1[8] 35.55
## m1[9] 39.48
## m1[10] 43.41
## y1[1] 3.90
## y1[2] -3.07
## y1[3] 23.56
## y1[4] 20.25
## y1[5] 13.94
## y1[6] 31.82
## y1[7] 27.56
## y1[8] 39.82
## y1[9] 45.52
## y1[10] 39.72
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -41.80 -41.43 1.29 1.08 -44.34 -40.36 1.00 464 605
## b0 7.01 6.84 2.93 2.80 2.46 12.07 1.01 297 327
## b1 1.85 1.86 0.25 0.23 1.43 2.24 1.01 348 406
## s 5.59 5.45 1.03 0.97 4.17 7.47 1.00 1031 1048
## m0[1] 28.99 28.99 1.29 1.28 26.85 31.12 1.00 1605 1744
## m0[2] 20.94 20.91 1.50 1.42 18.53 23.52 1.01 396 511
## m0[3] 38.68 38.75 2.00 1.88 35.29 41.75 1.00 990 1113
## m0[4] 15.83 15.78 1.95 1.85 12.77 19.14 1.01 320 363
## m0[5] 15.94 15.89 1.94 1.83 12.90 19.22 1.01 321 367
## m0[6] 16.73 16.68 1.86 1.74 13.81 19.92 1.01 327 356
## m0[7] 37.50 37.59 1.88 1.79 34.32 40.40 1.00 1109 1349
## m0[8] 38.04 38.12 1.93 1.82 34.76 41.01 1.00 1052 1235
## m0[9] 24.02 24.02 1.33 1.26 21.87 26.27 1.00 555 609
## m0[10] 21.79 21.78 1.44 1.38 19.40 24.23 1.01 425 557
## m0[11] 15.85 15.80 1.95 1.85 12.79 19.15 1.01 320 363
## m0[12] 30.29 30.31 1.34 1.32 28.07 32.46 1.00 1966 1718
## m0[13] 27.12 27.11 1.27 1.25 24.99 29.22 1.00 1016 1329
## m0[14] 32.99 33.04 1.49 1.45 30.52 35.31 1.00 2026 1324
## m0[15] 39.31 39.40 2.06 1.93 35.75 42.48 1.00 935 967
## m0[16] 8.02 7.87 2.81 2.69 3.66 12.95 1.01 297 327
## m0[17] 31.08 31.11 1.38 1.35 28.77 33.28 1.00 2129 1643
## m0[18] 28.43 28.42 1.28 1.26 26.31 30.54 1.00 1426 1750
## m0[19] 30.52 30.54 1.35 1.33 28.27 32.68 1.00 2017 1718
## m0[20] 43.38 43.45 2.51 2.34 39.10 47.25 1.00 703 744
## y0[1] 29.02 28.98 5.79 5.60 19.65 38.30 1.00 2187 2011
## y0[2] 20.85 20.80 5.80 5.47 11.59 30.23 1.00 1631 1518
## y0[3] 38.75 38.67 6.08 5.84 28.99 48.64 1.00 1714 1709
## y0[4] 15.73 15.75 5.95 5.73 6.16 25.49 1.00 1207 1483
## y0[5] 16.05 16.11 6.03 5.85 6.26 26.11 1.00 1316 1420
## y0[6] 16.62 16.76 5.95 5.84 6.75 26.05 1.00 1210 1243
## y0[7] 37.60 37.52 6.06 5.85 27.80 47.60 1.00 1974 1988
## y0[8] 38.21 38.36 5.93 5.70 28.29 47.64 1.00 1978 2007
## y0[9] 24.14 24.12 6.06 5.61 14.42 34.05 1.00 1642 1881
## y0[10] 21.71 21.82 5.72 5.50 12.67 30.95 1.00 1742 1671
## y0[11] 15.77 15.82 6.01 5.76 5.99 25.68 1.00 1272 1745
## y0[12] 30.18 30.06 5.80 5.46 20.70 39.74 1.00 2002 1597
## y0[13] 27.09 27.05 5.90 5.53 17.38 36.61 1.00 2160 2054
## y0[14] 32.89 32.97 5.85 5.70 23.30 42.76 1.00 1965 1900
## y0[15] 39.41 39.37 6.13 5.79 28.97 49.30 1.00 1998 1834
## y0[16] 8.04 7.96 6.24 6.05 -2.15 18.16 1.01 951 1250
## y0[17] 31.15 31.10 5.88 5.65 21.69 40.58 1.00 1905 1931
## y0[18] 28.35 28.35 5.82 5.68 18.76 37.99 1.00 1965 1894
## y0[19] 30.63 30.55 5.68 5.43 21.39 39.89 1.00 2261 1797
## y0[20] 43.35 43.43 6.25 5.84 33.10 53.53 1.00 1570 1729
## m1[1] 8.02 7.87 2.81 2.69 3.66 12.95 1.01 297 327
## m1[2] 11.95 11.87 2.36 2.28 8.16 16.05 1.01 302 336
## m1[3] 15.88 15.82 1.94 1.84 12.82 19.17 1.01 320 363
## m1[4] 19.81 19.78 1.59 1.49 17.30 22.55 1.01 370 437
## m1[5] 23.73 23.73 1.34 1.28 21.55 26.00 1.00 532 620
## m1[6] 27.66 27.65 1.27 1.26 25.54 29.75 1.00 1163 1473
## m1[7] 31.59 31.64 1.40 1.38 29.22 33.85 1.00 2171 1588
## m1[8] 35.52 35.61 1.70 1.62 32.70 38.17 1.00 1400 1363
## m1[9] 39.45 39.54 2.08 1.95 35.86 42.65 1.00 924 973
## m1[10] 43.38 43.45 2.51 2.34 39.10 47.25 1.00 703 744
## y1[1] 8.01 7.95 6.40 6.12 -2.19 18.58 1.00 1112 1455
## y1[2] 11.78 11.84 6.00 5.66 2.02 21.46 1.00 1420 1515
## y1[3] 15.73 15.84 5.99 5.96 6.01 25.56 1.00 1417 1589
## y1[4] 19.70 19.58 6.00 5.86 10.10 29.44 1.00 1540 2003
## y1[5] 23.88 23.87 5.84 5.80 14.27 33.16 1.00 1892 1930
## y1[6] 27.66 27.69 5.88 5.45 18.33 37.77 1.00 2068 1866
## y1[7] 31.45 31.72 5.85 5.79 21.69 40.75 1.00 2001 1971
## y1[8] 35.65 35.40 5.85 5.60 26.09 45.50 1.00 1879 2014
## y1[9] 39.36 39.43 6.06 5.74 29.60 49.02 1.00 1815 1910
## y1[10] 43.14 43.12 6.12 5.87 33.37 53.30 1.00 1480 1703
quadratic regression
y=b0+b2(x-b1)**2 ex8-1.stan
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real b2;
real<lower=0> s;
}
model{
y~normal(b0+b2*(x-b1)^2,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b2*(x[i]-b1)^2;
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b2*(x1[i]-b1)^2;
y1[i]=normal_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,0.5*(x-4)**2+5,1)
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-1.stan')
fn(mdl,data)## Initial log joint probability = -137.406
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 37 -9.26199 0.0001523 0.00171394 1 1 47
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -9.26
## b0 4.16
## b1 3.99
## b2 0.54
## s 0.96
## m0[1] 5.55
## m0[2] 5.28
## m0[3] 15.12
## m0[4] 14.50
## m0[5] 9.77
## m0[6] 11.75
## m0[7] 5.19
## m0[8] 10.77
## m0[9] 5.44
## m0[10] 7.54
## m0[11] 17.58
## m0[12] 8.87
## m0[13] 6.28
## m0[14] 4.58
## m0[15] 6.33
## m0[16] 8.61
## m0[17] 6.28
## m0[18] 9.71
## m0[19] 10.38
## m0[20] 14.79
## y0[1] 7.37
## y0[2] 6.68
## y0[3] 15.33
## y0[4] 13.88
## y0[5] 9.92
## y0[6] 12.58
## y0[7] 4.18
## y0[8] 11.46
## y0[9] 3.99
## y0[10] 9.51
## y0[11] 18.57
## y0[12] 8.92
## y0[13] 5.52
## y0[14] 5.61
## y0[15] 6.28
## y0[16] 8.05
## y0[17] 7.76
## y0[18] 9.14
## y0[19] 11.20
## y0[20] 13.57
## m1[1] 11.75
## m1[2] 8.33
## m1[3] 5.93
## m1[4] 4.54
## m1[5] 4.17
## m1[6] 4.82
## m1[7] 6.48
## m1[8] 9.16
## m1[9] 12.86
## m1[10] 17.58
## y1[1] 12.30
## y1[2] 8.88
## y1[3] 3.58
## y1[4] 3.38
## y1[5] 3.67
## y1[6] 4.19
## y1[7] 4.92
## y1[8] 8.97
## y1[9] 13.30
## y1[10] 17.73
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -11.46 -11.11 1.55 1.32 -14.37 -9.65 1.00 546 851
## b0 4.15 4.16 0.42 0.41 3.42 4.82 1.01 375 405
## b1 3.98 3.98 0.09 0.08 3.84 4.12 1.00 1580 1017
## b2 0.54 0.54 0.04 0.04 0.48 0.61 1.01 517 620
## s 1.14 1.10 0.22 0.20 0.85 1.54 1.01 879 889
## m0[1] 5.55 5.55 0.39 0.38 4.91 6.17 1.01 410 455
## m0[2] 5.27 5.28 0.37 0.37 4.65 5.87 1.01 443 533
## m0[3] 15.11 15.11 0.50 0.48 14.29 15.92 1.00 1567 1212
## m0[4] 14.49 14.49 0.47 0.45 13.73 15.25 1.00 1780 1270
## m0[5] 9.76 9.75 0.47 0.45 9.02 10.52 1.00 1617 1336
## m0[6] 11.74 11.74 0.59 0.56 10.83 12.69 1.00 1599 1337
## m0[7] 5.19 5.19 0.40 0.39 4.52 5.83 1.01 398 464
## m0[8] 10.76 10.77 0.34 0.33 10.22 11.31 1.00 1819 1658
## m0[9] 5.43 5.44 0.37 0.37 4.82 6.02 1.01 458 527
## m0[10] 7.53 7.53 0.34 0.33 6.98 8.08 1.01 550 767
## m0[11] 17.56 17.57 0.64 0.59 16.51 18.58 1.00 1123 1086
## m0[12] 8.86 8.85 0.42 0.40 8.19 9.58 1.00 1442 1319
## m0[13] 6.28 6.28 0.37 0.36 5.69 6.88 1.01 442 555
## m0[14] 4.58 4.60 0.39 0.39 3.91 5.21 1.01 394 448
## m0[15] 6.33 6.32 0.37 0.36 5.74 6.93 1.01 445 575
## m0[16] 8.60 8.60 0.32 0.33 8.07 9.13 1.01 744 1020
## m0[17] 6.27 6.27 0.37 0.36 5.68 6.87 1.01 442 555
## m0[18] 9.70 9.70 0.46 0.44 8.97 10.46 1.00 1611 1336
## m0[19] 10.37 10.36 0.50 0.48 9.58 11.19 1.00 1661 1267
## m0[20] 14.78 14.78 0.49 0.46 13.98 15.56 1.00 1676 1244
## y0[1] 5.59 5.63 1.22 1.16 3.62 7.58 1.00 1337 1911
## y0[2] 5.27 5.28 1.22 1.13 3.28 7.25 1.00 1542 1712
## y0[3] 15.12 15.14 1.25 1.19 13.07 17.15 1.00 1980 2017
## y0[4] 14.46 14.49 1.24 1.20 12.32 16.40 1.00 2054 1633
## y0[5] 9.78 9.76 1.28 1.21 7.74 11.89 1.00 1927 1834
## y0[6] 11.73 11.74 1.26 1.23 9.72 13.74 1.00 1944 1758
## y0[7] 5.19 5.19 1.25 1.16 3.17 7.16 1.00 1548 1386
## y0[8] 10.82 10.84 1.21 1.13 8.87 12.79 1.00 2122 1913
## y0[9] 5.40 5.45 1.22 1.17 3.35 7.37 1.00 1556 1688
## y0[10] 7.54 7.49 1.21 1.14 5.61 9.60 1.00 1450 1740
## y0[11] 17.56 17.59 1.30 1.21 15.37 19.62 1.00 1928 1713
## y0[12] 8.87 8.90 1.25 1.17 6.74 10.89 1.00 1779 1756
## y0[13] 6.30 6.33 1.18 1.12 4.36 8.21 1.00 1982 1880
## y0[14] 4.59 4.61 1.21 1.12 2.65 6.60 1.00 1457 1684
## y0[15] 6.32 6.31 1.19 1.14 4.42 8.26 1.00 1238 1745
## y0[16] 8.58 8.57 1.20 1.12 6.72 10.51 1.00 1713 1685
## y0[17] 6.27 6.30 1.23 1.12 4.30 8.29 1.00 1707 1758
## y0[18] 9.69 9.69 1.21 1.15 7.73 11.65 1.00 1937 1664
## y0[19] 10.33 10.34 1.25 1.15 8.26 12.38 1.00 1825 1750
## y0[20] 14.74 14.76 1.23 1.16 12.73 16.75 1.00 1925 1882
## m1[1] 11.74 11.74 0.59 0.56 10.83 12.69 1.00 1599 1337
## m1[2] 8.32 8.32 0.40 0.38 7.69 9.01 1.00 1248 1249
## m1[3] 5.92 5.93 0.36 0.36 5.33 6.50 1.01 520 573
## m1[4] 4.54 4.55 0.40 0.40 3.87 5.18 1.01 392 443
## m1[5] 4.17 4.17 0.42 0.42 3.44 4.84 1.01 375 419
## m1[6] 4.81 4.82 0.41 0.41 4.12 5.47 1.01 388 431
## m1[7] 6.48 6.48 0.36 0.35 5.89 7.07 1.01 454 606
## m1[8] 9.16 9.16 0.32 0.32 8.63 9.68 1.00 910 1279
## m1[9] 12.85 12.85 0.40 0.38 12.20 13.50 1.00 2073 1664
## m1[10] 17.56 17.57 0.64 0.59 16.51 18.58 1.00 1123 1086
## y1[1] 11.74 11.76 1.30 1.27 9.61 13.81 1.00 2003 1870
## y1[2] 8.28 8.27 1.21 1.12 6.30 10.28 1.00 1740 1876
## y1[3] 5.92 5.88 1.24 1.16 4.03 8.00 1.00 1545 1858
## y1[4] 4.56 4.54 1.22 1.16 2.57 6.62 1.00 1419 1688
## y1[5] 4.16 4.16 1.22 1.20 2.18 6.12 1.00 1443 1649
## y1[6] 4.83 4.81 1.23 1.21 2.76 6.87 1.00 1413 1642
## y1[7] 6.48 6.53 1.21 1.17 4.45 8.47 1.00 1453 1719
## y1[8] 9.15 9.15 1.22 1.19 7.12 11.18 1.00 1816 1813
## y1[9] 12.85 12.82 1.21 1.17 10.89 14.82 1.00 1814 1855
## y1[10] 17.56 17.56 1.33 1.31 15.39 19.67 1.00 1840 1578
both log regression
log y=b0+b1*log x x,y>0
y=exp b0 * x**b1ex8-2.stan
data{
int N;
vector<lower=0>[N] x;
vector<lower=0>[N] y;
int N1;
vector<lower=0>[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~lognormal(b0+b1*log(x),s);
}
generated quantities{
vector[N] m0;
vector<lower=0>[N] y0;
for(i in 1:N){
m0[i]=b0+b1*log(x[i]);
y0[i]=lognormal_rng(m0[i],s);
}
vector[N1] m1;
vector<lower=0>[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*log(x1[i]);
y1[i]=lognormal_rng(m1[i],s);
}
}n=20
x=runif(n,0,10)
y=exp(rnorm(n,log(x)*2+1,1))
grid.arrange(qplot(x,y),
qplot(log(x),log(y)),
ncol=2)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -5181.34
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 15 -29.9396 0.000262506 0.000132818 1 1 30
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -29.94
## b0 1.05
## b1 2.16
## s 1.08
## m0[1] 5.31
## m0[2] 5.60
## m0[3] 1.36
## m0[4] 5.97
## m0[5] 2.09
## m0[6] 5.86
## m0[7] 2.52
## m0[8] 5.59
## m0[9] 6.01
## m0[10] 3.37
## m0[11] 4.61
## m0[12] 4.58
## m0[13] 5.06
## m0[14] 5.62
## m0[15] 5.45
## m0[16] 5.27
## m0[17] -2.68
## m0[18] 1.15
## m0[19] 4.28
## m0[20] 0.65
## y0[1] 40.62
## y0[2] 343.62
## y0[3] 1.00
## y0[4] 48.48
## y0[5] 3.30
## y0[6] 59.33
## y0[7] 83.10
## y0[8] 1195.92
## y0[9] 1455.42
## y0[10] 23.08
## y0[11] 163.65
## y0[12] 525.74
## y0[13] 90.22
## y0[14] 135.41
## y0[15] 276.20
## y0[16] 101.27
## y0[17] 0.14
## y0[18] 3.47
## y0[19] 199.75
## y0[20] 2.77
## m1[1] -2.68
## m1[2] 1.55
## m1[3] 2.89
## m1[4] 3.71
## m1[5] 4.30
## m1[6] 4.77
## m1[7] 5.15
## m1[8] 5.48
## m1[9] 5.76
## m1[10] 6.01
## y1[1] 0.03
## y1[2] 7.12
## y1[3] 21.24
## y1[4] 60.92
## y1[5] 78.03
## y1[6] 412.04
## y1[7] 70.63
## y1[8] 2071.88
## y1[9] 124.94
## y1[10] 250.87mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.mcmc$metadata()$stan_variables## [1] "lp__" "b0" "b1" "s" "m0" "y0" "m1" "y1"seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -31.50 -31.16 1.34 1.08 -34.22 -30.05 1.00 706 1024
## b0 1.04 1.06 0.47 0.46 0.25 1.76 1.01 393 686
## b1 2.17 2.16 0.28 0.27 1.73 2.62 1.01 478 861
## s 1.23 1.20 0.23 0.21 0.93 1.65 1.00 1385 1425
## m0[1] 5.31 5.31 0.32 0.31 4.78 5.82 1.01 1952 1457
## m0[2] 5.60 5.60 0.34 0.33 5.03 6.15 1.00 1833 1394
## m0[3] 1.35 1.37 0.43 0.42 0.61 2.03 1.01 397 713
## m0[4] 5.98 5.98 0.37 0.37 5.36 6.57 1.00 1569 1457
## m0[5] 2.08 2.10 0.37 0.36 1.48 2.67 1.01 412 763
## m0[6] 5.87 5.87 0.36 0.35 5.26 6.44 1.00 1644 1374
## m0[7] 2.51 2.52 0.33 0.33 1.98 3.05 1.01 431 855
## m0[8] 5.60 5.60 0.34 0.33 5.03 6.14 1.00 1836 1394
## m0[9] 6.02 6.02 0.37 0.37 5.40 6.61 1.00 1538 1435
## m0[10] 3.37 3.37 0.28 0.27 2.91 3.83 1.01 575 1028
## m0[11] 4.61 4.61 0.28 0.27 4.15 5.06 1.00 1448 1665
## m0[12] 4.58 4.58 0.28 0.27 4.13 5.03 1.00 1416 1690
## m0[13] 5.07 5.07 0.30 0.29 4.56 5.56 1.00 1888 1586
## m0[14] 5.63 5.63 0.34 0.33 5.05 6.18 1.00 1817 1394
## m0[15] 5.46 5.46 0.33 0.32 4.91 5.99 1.00 1917 1437
## m0[16] 5.27 5.27 0.31 0.30 4.75 5.78 1.01 1969 1459
## m0[17] -2.72 -2.68 0.91 0.89 -4.24 -1.34 1.01 403 721
## m0[18] 1.14 1.16 0.46 0.45 0.37 1.85 1.01 394 725
## m0[19] 4.28 4.28 0.27 0.26 3.83 4.71 1.01 1067 1356
## m0[20] 0.64 0.66 0.51 0.50 -0.24 1.41 1.01 390 671
## y0[1] 483.09 201.70 1255.66 207.38 25.15 1601.75 1.00 1877 1932
## y0[2] 648.86 274.92 2153.91 290.67 33.93 2209.18 1.00 2138 2044
## y0[3] 9.08 3.86 20.12 4.02 0.44 31.41 1.00 1283 1966
## y0[4] 974.43 382.68 2611.32 402.41 44.02 3163.10 1.00 2038 1846
## y0[5] 18.92 8.32 65.03 8.79 0.96 61.62 1.00 1619 1973
## y0[6] 832.19 359.83 1786.01 375.18 38.78 2957.74 1.00 1934 1712
## y0[7] 30.65 12.70 70.65 13.96 1.40 105.57 1.00 1624 1764
## y0[8] 806.97 260.03 9360.08 270.99 29.11 2067.98 1.00 2048 1696
## y0[9] 1063.43 434.67 2777.23 450.57 56.50 3707.01 1.00 1868 1773
## y0[10] 68.00 29.38 144.46 30.71 3.52 242.03 1.00 1677 1970
## y0[11] 238.16 99.98 731.96 103.69 12.46 757.54 1.00 1978 1856
## y0[12] 229.47 94.96 482.03 95.30 12.06 821.65 1.00 2101 1799
## y0[13] 366.09 160.51 806.35 168.53 18.40 1382.51 1.00 2035 1945
## y0[14] 702.27 269.73 2291.80 278.49 32.49 2346.89 1.00 1837 1647
## y0[15] 531.43 222.82 1377.34 233.42 25.85 1760.23 1.00 2065 1891
## y0[16] 499.38 192.34 1380.56 198.56 24.29 1770.73 1.00 2169 1927
## y0[17] 0.24 0.07 0.97 0.08 0.01 0.85 1.01 814 1322
## y0[18] 7.58 3.25 17.59 3.53 0.35 27.53 1.00 1473 1831
## y0[19] 160.03 70.96 283.27 74.62 7.78 571.37 1.00 1623 1875
## y0[20] 4.50 1.84 10.33 2.05 0.19 15.70 1.00 1720 1593
## m1[1] -2.72 -2.68 0.91 0.89 -4.24 -1.34 1.01 403 721
## m1[2] 1.54 1.55 0.42 0.41 0.84 2.19 1.01 401 708
## m1[3] 2.88 2.89 0.30 0.30 2.39 3.39 1.01 467 992
## m1[4] 3.71 3.71 0.27 0.26 3.26 4.15 1.01 702 1101
## m1[5] 4.30 4.30 0.27 0.26 3.86 4.74 1.01 1088 1425
## m1[6] 4.77 4.77 0.29 0.28 4.30 5.24 1.00 1621 1619
## m1[7] 5.16 5.16 0.31 0.29 4.64 5.66 1.00 1930 1587
## m1[8] 5.48 5.48 0.33 0.32 4.93 6.02 1.00 1905 1437
## m1[9] 5.77 5.77 0.35 0.35 5.17 6.33 1.00 1718 1375
## m1[10] 6.02 6.02 0.37 0.37 5.40 6.61 1.00 1538 1435
## y1[1] 0.23 0.06 0.75 0.08 0.01 0.75 1.01 977 1254
## y1[2] 14.27 4.60 136.39 4.64 0.55 40.32 1.00 1531 1935
## y1[3] 46.27 17.50 125.87 18.73 2.01 162.26 1.00 1652 1892
## y1[4] 97.46 39.66 238.91 41.78 4.62 305.60 1.00 2085 1918
## y1[5] 176.58 77.52 418.60 81.21 8.36 596.34 1.00 1951 1972
## y1[6] 307.66 110.52 982.07 117.52 13.80 968.40 1.00 2087 1988
## y1[7] 412.31 157.17 1033.51 165.26 20.78 1434.79 1.00 1995 1845
## y1[8] 600.40 253.31 1467.35 268.70 32.10 2043.06 1.00 1944 1890
## y1[9] 805.42 322.38 1953.98 325.43 41.24 2824.73 1.00 1996 1914
## y1[10] 986.18 394.15 2306.50 414.51 50.62 3518.26 1.00 2016 1848
y0=mcmc$draws('y0')
smy0=summary(y0)
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(log(y)-log(smy0$median),xlab='obs.-prd.',main='residual')
density(log(y)-log(smy0$median)) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=exp(m)),xy,col='black')
qplot(log(x),log(y),col=I('red'))+
geom_line(aes(x=log(x),y=ml5),xy,col='darkgray')+
geom_line(aes(x=log(x),y=mu5),xy,col='darkgray')+
geom_line(aes(x=log(x),y=log(yl5)),xy,col='lightgray')+
geom_line(aes(x=log(x),y=log(yu5)),xy,col='lightgray')+
geom_line(aes(x=log(x),y=m),xy,col='black')
ex8-3
exponential increasing/decreasing
y=b0* exp b1* x -> y~N(b0* exp(b1*x),s)
log y=log b0+b1* x -> y~logN(log b0 +b1*x,s)
x,y>0,b0>0
(x=0,y=b0)
b1>0 x->Infinity,y->Infinity
b1<0 x->Infinity,y->+0 n=20
x=runif(n,0,5)
y=rnorm(n,10*exp(-2*x),0.5)
grid.arrange(qplot(x,y),
qplot(x,log(y)),
ncol=2)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)y=b0* exp b1* x -> y~N(b0* exp(b1*x),s) ex8-3-1.stan
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0> b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0*exp(b1*x),s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0*exp(b1*x[i]);
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0*exp(b1*x1[i]);
y1[i]=normal_rng(m1[i],s);
}
}mdl=cmdstan_model('./ex8-3-1.stan')
fn(mdl,data)## Initial log joint probability = -24.6144
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 19 10.8609 4.19475e-05 0.00133438 1 1 30
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ 10.86
## b0 15.08
## b1 -2.34
## s 0.35
## m0[1] 0.00
## m0[2] 0.02
## m0[3] 3.50
## m0[4] 0.00
## m0[5] 0.00
## m0[6] 0.03
## m0[7] 0.00
## m0[8] 0.37
## m0[9] 0.00
## m0[10] 1.08
## m0[11] 0.11
## m0[12] 0.00
## m0[13] 0.12
## m0[14] 1.54
## m0[15] 0.69
## m0[16] 1.89
## m0[17] 1.92
## m0[18] 0.00
## m0[19] 0.68
## m0[20] 0.09
## y0[1] -0.44
## y0[2] 0.67
## y0[3] 3.91
## y0[4] 0.09
## y0[5] 0.01
## y0[6] 0.09
## y0[7] 0.07
## y0[8] 0.49
## y0[9] -0.04
## y0[10] 0.58
## y0[11] 0.22
## y0[12] -0.15
## y0[13] 0.18
## y0[14] 1.26
## y0[15] 1.58
## y0[16] 1.95
## y0[17] 1.97
## y0[18] 0.13
## y0[19] 0.96
## y0[20] -0.03
## m1[1] 3.50
## m1[2] 1.28
## m1[3] 0.47
## m1[4] 0.17
## m1[5] 0.06
## m1[6] 0.02
## m1[7] 0.01
## m1[8] 0.00
## m1[9] 0.00
## m1[10] 0.00
## y1[1] 3.74
## y1[2] 1.40
## y1[3] 0.11
## y1[4] 0.13
## y1[5] 0.17
## y1[6] -0.01
## y1[7] 0.07
## y1[8] 0.16
## y1[9] -0.64
## y1[10] -0.03
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 0.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 12.42 12.77 1.40 1.13 9.79 13.93 1.00 477 436
## b0 19.45 17.18 9.28 5.97 10.39 35.30 1.02 287 234
## b1 -2.58 -2.49 0.50 0.44 -3.47 -1.90 1.01 285 255
## s 0.41 0.39 0.08 0.07 0.30 0.55 1.00 604 731
## m0[1] 0.00 0.00 0.00 0.00 0.00 0.01 1.01 289 273
## m0[2] 0.02 0.02 0.02 0.01 0.00 0.06 1.01 291 302
## m0[3] 3.60 3.58 0.40 0.39 2.97 4.27 1.00 489 459
## m0[4] 0.00 0.00 0.00 0.00 0.00 0.01 1.01 289 285
## m0[5] 0.00 0.00 0.00 0.00 0.00 0.01 1.01 289 285
## m0[6] 0.03 0.02 0.02 0.02 0.00 0.07 1.01 291 302
## m0[7] 0.00 0.00 0.00 0.00 0.00 0.00 1.01 288 274
## m0[8] 0.32 0.32 0.12 0.12 0.14 0.53 1.01 305 325
## m0[9] 0.00 0.00 0.00 0.00 0.00 0.00 1.01 288 274
## m0[10] 1.00 1.01 0.19 0.19 0.68 1.30 1.01 360 335
## m0[11] 0.10 0.09 0.06 0.05 0.02 0.20 1.01 295 294
## m0[12] 0.00 0.00 0.00 0.00 0.00 0.00 1.01 288 274
## m0[13] 0.10 0.09 0.06 0.05 0.03 0.21 1.01 295 294
## m0[14] 1.46 1.47 0.20 0.19 1.13 1.78 1.01 472 388
## m0[15] 0.62 0.62 0.17 0.16 0.35 0.89 1.01 321 339
## m0[16] 1.83 1.83 0.20 0.18 1.49 2.14 1.00 694 598
## m0[17] 1.85 1.86 0.20 0.18 1.52 2.17 1.00 722 617
## m0[18] 0.00 0.00 0.00 0.00 0.00 0.00 1.01 288 274
## m0[19] 0.61 0.61 0.16 0.16 0.34 0.88 1.01 321 339
## m0[20] 0.08 0.07 0.05 0.04 0.02 0.17 1.01 294 292
## y0[1] 0.00 0.00 0.40 0.39 -0.66 0.65 1.00 1735 1797
## y0[2] 0.00 -0.01 0.43 0.41 -0.69 0.70 1.00 1870 1893
## y0[3] 3.61 3.60 0.58 0.55 2.70 4.54 1.00 895 972
## y0[4] -0.01 0.00 0.42 0.40 -0.68 0.66 1.00 2075 1995
## y0[5] -0.01 -0.01 0.41 0.39 -0.68 0.66 1.00 1998 2012
## y0[6] 0.04 0.04 0.41 0.38 -0.62 0.72 1.00 1920 2053
## y0[7] -0.01 -0.01 0.41 0.39 -0.69 0.66 1.00 2064 2052
## y0[8] 0.33 0.34 0.41 0.40 -0.36 0.96 1.00 1646 1809
## y0[9] -0.01 0.00 0.41 0.38 -0.67 0.66 1.00 1992 1591
## y0[10] 0.97 0.98 0.45 0.42 0.22 1.72 1.00 1148 1248
## y0[11] 0.10 0.11 0.41 0.39 -0.59 0.79 1.00 1886 1887
## y0[12] 0.00 0.00 0.42 0.39 -0.68 0.68 1.00 2052 1940
## y0[13] 0.12 0.12 0.42 0.40 -0.56 0.83 1.00 1918 1766
## y0[14] 1.46 1.46 0.47 0.43 0.67 2.22 1.00 1664 1895
## y0[15] 0.61 0.61 0.43 0.42 -0.09 1.32 1.00 1410 1696
## y0[16] 1.82 1.84 0.47 0.44 1.03 2.57 1.00 1433 1785
## y0[17] 1.85 1.86 0.45 0.45 1.08 2.57 1.00 1784 1512
## y0[18] 0.01 0.01 0.42 0.40 -0.67 0.72 1.00 2026 1893
## y0[19] 0.62 0.62 0.45 0.43 -0.13 1.33 1.00 1456 1280
## y0[20] 0.08 0.08 0.41 0.40 -0.59 0.75 1.00 2209 1837
## m1[1] 3.60 3.58 0.40 0.39 2.97 4.27 1.00 489 459
## m1[2] 1.20 1.20 0.19 0.19 0.87 1.50 1.01 394 348
## m1[3] 0.42 0.41 0.14 0.14 0.20 0.65 1.01 309 340
## m1[4] 0.15 0.14 0.08 0.07 0.05 0.29 1.01 297 303
## m1[5] 0.06 0.05 0.04 0.03 0.01 0.13 1.01 293 292
## m1[6] 0.02 0.02 0.02 0.01 0.00 0.06 1.01 291 302
## m1[7] 0.01 0.01 0.01 0.01 0.00 0.02 1.01 290 285
## m1[8] 0.00 0.00 0.01 0.00 0.00 0.01 1.01 289 285
## m1[9] 0.00 0.00 0.00 0.00 0.00 0.00 1.01 288 274
## m1[10] 0.00 0.00 0.00 0.00 0.00 0.00 1.01 288 274
## y1[1] 3.59 3.60 0.57 0.57 2.67 4.53 1.00 864 908
## y1[2] 1.20 1.21 0.47 0.45 0.43 1.94 1.00 1374 1336
## y1[3] 0.41 0.42 0.44 0.41 -0.33 1.13 1.00 1519 1678
## y1[4] 0.15 0.15 0.42 0.39 -0.52 0.83 1.00 1974 2055
## y1[5] 0.06 0.05 0.42 0.41 -0.64 0.74 1.00 1975 1562
## y1[6] 0.03 0.04 0.42 0.42 -0.68 0.69 1.00 1945 1811
## y1[7] 0.01 0.00 0.42 0.40 -0.69 0.70 1.00 2025 1625
## y1[8] 0.03 0.04 0.42 0.40 -0.65 0.72 1.00 1938 1840
## y1[9] 0.02 0.02 0.41 0.40 -0.61 0.72 1.00 2338 1857
## y1[10] 0.00 0.02 0.41 0.40 -0.68 0.64 1.00 1934 1895
log y=log b0+b1* x -> y~logN(log b0 +b1*x,s)ex8-3-2.stan
data{
int N;
vector<lower=0>[N] x;
vector<lower=0>[N] y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0> b0;
real<lower=-10,upper=10> b1;
real<lower=0> s;
}
model{
y~lognormal(log(b0)+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=log(b0)+b1*x[i];
y0[i]=lognormal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=log(b0)+b1*x1[i];
y1[i]=lognormal_rng(m1[i],s);
}
}n=20
x=runif(n,0,5)
y=rlnorm(n,log(10)-2*x,0.5)
grid.arrange(qplot(x,y),
qplot(x,log(y)),
ncol=2)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-3-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -14383.4
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 12 -79.5531 0.00216098 0.0003441 0.9588 0.9588 31
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -79.55
## b0 17121900000.00
## b1 -10.00
## s 12.92
## m0[1] 16.93
## m0[2] -25.28
## m0[3] 10.32
## m0[4] -17.48
## m0[5] -25.54
## m0[6] 8.44
## m0[7] 23.15
## m0[8] -18.20
## m0[9] -21.39
## m0[10] -16.04
## m0[11] -11.34
## m0[12] 9.98
## m0[13] -18.46
## m0[14] 6.64
## m0[15] -20.83
## m0[16] 1.48
## m0[17] -10.69
## m0[18] 18.36
## m0[19] 18.71
## m0[20] 4.04
## y0[1] 1.37
## y0[2] 2.88
## y0[3] 168823000000.00
## y0[4] 0.00
## y0[5] 0.00
## y0[6] 6574500.00
## y0[7] 0.00
## y0[8] 0.00
## y0[9] 0.00
## y0[10] 0.00
## y0[11] 0.02
## y0[12] 30023400.00
## y0[13] 0.09
## y0[14] 286488.00
## y0[15] 0.00
## y0[16] 8.80
## y0[17] 0.00
## y0[18] 14414600000.00
## y0[19] 13937800000.00
## y0[20] 0.00
## m1[1] 23.15
## m1[2] 17.74
## m1[3] 12.33
## m1[4] 6.92
## m1[5] 1.51
## m1[6] -3.90
## m1[7] -9.31
## m1[8] -14.72
## m1[9] -20.13
## m1[10] -25.54
## y1[1] 67628300.00
## y1[2] 13202.10
## y1[3] 792279000.00
## y1[4] 10221900000000000000.00
## y1[5] 0.00
## y1[6] 0.03
## y1[7] 0.00
## y1[8] 2.50
## y1[9] 0.00
## y1[10] 0.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.mcmc$metadata()$stan_variables## [1] "lp__" "b0" "b1" "s" "m0" "y0" "m1" "y1"seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -7.40 -7.08 1.26 1.04 -9.78 -6.05 1.01 654 829
## b0 9.36 9.09 1.72 1.55 6.96 12.49 1.01 650 473
## b1 -2.07 -2.07 0.06 0.05 -2.17 -1.98 1.00 831 700
## s 0.41 0.40 0.07 0.06 0.31 0.55 1.01 686 664
## m0[1] 0.85 0.84 0.15 0.14 0.62 1.10 1.01 655 463
## m0[2] -7.89 -7.88 0.16 0.15 -8.14 -7.63 1.00 1811 1327
## m0[3] -0.52 -0.53 0.12 0.12 -0.71 -0.31 1.00 699 525
## m0[4] -6.27 -6.27 0.13 0.12 -6.48 -6.07 1.00 2134 1590
## m0[5] -7.94 -7.94 0.16 0.16 -8.20 -7.68 1.00 1800 1326
## m0[6] -0.91 -0.91 0.11 0.11 -1.09 -0.71 1.00 726 597
## m0[7] 2.13 2.12 0.17 0.17 1.86 2.43 1.01 650 473
## m0[8] -6.42 -6.42 0.13 0.13 -6.63 -6.21 1.00 2111 1445
## m0[9] -7.08 -7.08 0.14 0.14 -7.31 -6.86 1.00 1980 1425
## m0[10] -5.97 -5.97 0.12 0.12 -6.17 -5.78 1.00 2165 1374
## m0[11] -5.00 -5.00 0.11 0.11 -5.17 -4.83 1.00 2153 1508
## m0[12] -0.59 -0.60 0.12 0.12 -0.78 -0.38 1.00 704 525
## m0[13] -6.47 -6.47 0.13 0.13 -6.69 -6.27 1.00 2102 1444
## m0[14] -1.28 -1.29 0.11 0.10 -1.45 -1.09 1.00 758 586
## m0[15] -6.96 -6.96 0.14 0.14 -7.19 -6.75 1.00 2004 1403
## m0[16] -2.35 -2.35 0.10 0.09 -2.51 -2.18 1.00 931 802
## m0[17] -4.87 -4.87 0.10 0.10 -5.04 -4.70 1.00 2117 1555
## m0[18] 1.14 1.14 0.15 0.15 0.90 1.41 1.01 652 475
## m0[19] 1.22 1.21 0.15 0.15 0.97 1.48 1.01 651 475
## m0[20] -1.82 -1.82 0.10 0.10 -1.98 -1.65 1.00 824 676
## y0[1] 2.62 2.34 1.24 0.96 1.19 5.03 1.00 1710 1618
## y0[2] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 2073 1550
## y0[3] 0.65 0.59 0.31 0.25 0.28 1.21 1.00 1690 1846
## y0[4] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 2103 1933
## y0[5] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 2050 1884
## y0[6] 0.45 0.40 0.21 0.16 0.20 0.85 1.00 1988 1970
## y0[7] 9.43 8.41 4.87 3.62 4.04 18.50 1.00 1352 1470
## y0[8] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1865 1955
## y0[9] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 2062 1945
## y0[10] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 1854 1714
## y0[11] 0.01 0.01 0.00 0.00 0.00 0.01 1.00 1960 1721
## y0[12] 0.60 0.55 0.27 0.22 0.28 1.07 1.00 1858 1855
## y0[13] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1962 1841
## y0[14] 0.30 0.27 0.13 0.11 0.14 0.54 1.00 1806 1893
## y0[15] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1992 1969
## y0[16] 0.11 0.10 0.05 0.04 0.05 0.19 1.00 1751 1839
## y0[17] 0.01 0.01 0.00 0.00 0.00 0.02 1.00 2095 2011
## y0[18] 3.49 3.13 1.79 1.32 1.42 6.65 1.00 1480 1608
## y0[19] 3.70 3.36 1.80 1.41 1.57 7.10 1.00 1420 1922
## y0[20] 0.18 0.16 0.08 0.06 0.08 0.33 1.00 1764 1567
## m1[1] 2.13 2.12 0.17 0.17 1.86 2.43 1.01 650 473
## m1[2] 1.02 1.01 0.15 0.15 0.78 1.27 1.01 653 463
## m1[3] -0.10 -0.11 0.13 0.12 -0.31 0.11 1.00 679 500
## m1[4] -1.22 -1.23 0.11 0.11 -1.39 -1.04 1.00 753 580
## m1[5] -2.34 -2.34 0.10 0.09 -2.50 -2.18 1.00 929 802
## m1[6] -3.46 -3.46 0.09 0.09 -3.62 -3.31 1.00 1398 1537
## m1[7] -4.58 -4.58 0.10 0.10 -4.74 -4.42 1.00 2020 1718
## m1[8] -5.70 -5.70 0.12 0.11 -5.89 -5.52 1.00 2187 1439
## m1[9] -6.82 -6.82 0.14 0.13 -7.04 -6.61 1.00 2035 1362
## m1[10] -7.94 -7.94 0.16 0.16 -8.20 -7.68 1.00 1800 1326
## y1[1] 9.21 8.34 4.59 3.53 3.88 16.99 1.00 1495 1939
## y1[2] 3.03 2.73 1.53 1.17 1.35 5.67 1.00 1864 1709
## y1[3] 0.99 0.90 0.47 0.38 0.43 1.82 1.00 1795 1820
## y1[4] 0.32 0.29 0.15 0.12 0.15 0.60 1.01 1791 1910
## y1[5] 0.11 0.10 0.05 0.04 0.05 0.20 1.00 1854 1970
## y1[6] 0.03 0.03 0.02 0.01 0.02 0.06 1.00 1597 1949
## y1[7] 0.01 0.01 0.01 0.00 0.01 0.02 1.00 1821 1651
## y1[8] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 1911 1919
## y1[9] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1757 2012
## y1[10] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1763 1889
y0=mcmc$draws('y0')
smy0=summary(y0)
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,log(y),col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=log(yl5)),xy,col='lightgray')+
geom_line(aes(x=x,y=log(yu5)),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
ex8-4
growth curve
y=b0* (1-exp(-b1* x)) -> y~N(1-exp(-b1*x),s)
x,y>0, b0,b1>0
(x=0,y=0), (x->Infinity,y->b0)data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0,upper=100> b0;
real<lower=0,upper=10> b1;
real<lower=0> s;
}
model{
y~normal(b0*(1-exp(-b1*x)),s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0*(1-exp(-b1*x[i]));
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0*(1-exp(-b1*x1[i]));
y1[i]=normal_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,10*(1-exp(-0.5*x)),1)
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-4.stan')
fn(mdl,data)## Initial log joint probability = -222558
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 -51.0573 0.000272808 0.000222331 0.9729 0.9729 23
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -51.06
## b0 0.00
## b1 0.00
## s 7.79
## m0[1] 0.00
## m0[2] 0.00
## m0[3] 0.00
## m0[4] 0.00
## m0[5] 0.00
## m0[6] 0.00
## m0[7] 0.00
## m0[8] 0.00
## m0[9] 0.00
## m0[10] 0.00
## m0[11] 0.00
## m0[12] 0.00
## m0[13] 0.00
## m0[14] 0.00
## m0[15] 0.00
## m0[16] 0.00
## m0[17] 0.00
## m0[18] 0.00
## m0[19] 0.00
## m0[20] 0.00
## y0[1] -13.21
## y0[2] -7.29
## y0[3] -8.49
## y0[4] -2.52
## y0[5] 14.14
## y0[6] 14.25
## y0[7] -10.18
## y0[8] -6.30
## y0[9] -1.78
## y0[10] 5.56
## y0[11] -9.93
## y0[12] -13.70
## y0[13] 3.12
## y0[14] 1.31
## y0[15] 4.51
## y0[16] 7.46
## y0[17] -5.70
## y0[18] -4.72
## y0[19] 3.04
## y0[20] 5.99
## m1[1] 0.00
## m1[2] 0.00
## m1[3] 0.00
## m1[4] 0.00
## m1[5] 0.00
## m1[6] 0.00
## m1[7] 0.00
## m1[8] 0.00
## m1[9] 0.00
## m1[10] 0.00
## y1[1] -10.61
## y1[2] 2.48
## y1[3] 0.98
## y1[4] -1.02
## y1[5] -6.68
## y1[6] -12.56
## y1[7] 11.68
## y1[8] 4.19
## y1[9] 3.50
## y1[10] 3.16
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -5.43 -5.09 1.31 1.07 -8.10 -4.01 1.01 688 795
## b0 10.53 10.52 0.47 0.43 9.82 11.31 1.00 923 1068
## b1 0.44 0.44 0.06 0.05 0.36 0.54 1.00 762 768
## s 0.89 0.87 0.16 0.14 0.67 1.17 1.01 570 499
## m0[1] 4.51 4.50 0.29 0.27 4.05 5.01 1.00 886 891
## m0[2] 10.27 10.27 0.36 0.36 9.69 10.85 1.00 1107 1174
## m0[3] 5.96 5.95 0.31 0.30 5.45 6.49 1.00 953 919
## m0[4] 1.84 1.83 0.16 0.14 1.60 2.10 1.00 841 863
## m0[5] 9.87 9.87 0.28 0.28 9.40 10.31 1.00 1689 1428
## m0[6] 8.17 8.17 0.26 0.26 7.73 8.60 1.00 1468 1153
## m0[7] 8.39 8.39 0.26 0.25 7.97 8.79 1.00 1628 1248
## m0[8] 4.72 4.71 0.30 0.28 4.25 5.23 1.00 892 891
## m0[9] 9.58 9.59 0.25 0.24 9.18 9.99 1.00 2130 1649
## m0[10] 10.27 10.27 0.36 0.36 9.68 10.85 1.00 1107 1174
## m0[11] 8.02 8.02 0.27 0.27 7.57 8.45 1.00 1375 1126
## m0[12] 1.29 1.28 0.11 0.10 1.12 1.48 1.00 836 840
## m0[13] 1.19 1.19 0.11 0.10 1.03 1.37 1.00 835 840
## m0[14] 10.23 10.23 0.35 0.35 9.66 10.79 1.00 1149 1189
## m0[15] 9.48 9.48 0.25 0.24 9.08 9.88 1.00 2259 1568
## m0[16] 9.96 9.97 0.29 0.29 9.48 10.43 1.00 1490 1451
## m0[17] 9.74 9.75 0.27 0.26 9.31 10.17 1.00 1886 1469
## m0[18] 8.82 8.82 0.24 0.23 8.42 9.21 1.00 2051 1648
## m0[19] 5.24 5.23 0.31 0.29 4.75 5.77 1.00 911 905
## m0[20] 4.42 4.42 0.29 0.27 3.97 4.92 1.00 883 891
## y0[1] 4.52 4.49 0.93 0.88 3.07 5.98 1.00 1661 1865
## y0[2] 10.28 10.27 0.96 0.94 8.72 11.86 1.00 1945 1711
## y0[3] 5.95 5.96 0.96 0.91 4.43 7.52 1.00 1807 1768
## y0[4] 1.85 1.83 0.92 0.89 0.38 3.36 1.00 1827 1802
## y0[5] 9.86 9.86 0.96 0.96 8.29 11.42 1.00 1888 1789
## y0[6] 8.15 8.15 0.92 0.85 6.62 9.65 1.00 1903 1931
## y0[7] 8.42 8.41 0.95 0.90 6.88 9.97 1.00 1956 1717
## y0[8] 4.71 4.72 0.95 0.89 3.15 6.27 1.00 1886 1648
## y0[9] 9.59 9.58 0.95 0.91 8.03 11.07 1.00 1986 1953
## y0[10] 10.27 10.26 0.98 0.94 8.68 11.88 1.00 1809 1711
## y0[11] 7.99 8.01 0.95 0.91 6.39 9.57 1.00 1859 1621
## y0[12] 1.30 1.29 0.92 0.86 -0.21 2.78 1.00 1827 1742
## y0[13] 1.21 1.25 0.91 0.91 -0.36 2.68 1.00 1967 1721
## y0[14] 10.23 10.21 0.99 0.96 8.64 11.86 1.00 1937 1736
## y0[15] 9.48 9.47 0.94 0.91 7.94 11.05 1.00 1949 1837
## y0[16] 9.99 9.95 0.93 0.90 8.47 11.51 1.00 1782 1729
## y0[17] 9.73 9.72 0.98 0.94 8.16 11.36 1.00 1926 1809
## y0[18] 8.81 8.79 0.95 0.91 7.25 10.35 1.00 1949 1708
## y0[19] 5.24 5.26 0.95 0.87 3.63 6.81 1.00 1780 1774
## y0[20] 4.38 4.39 0.97 0.94 2.83 5.91 1.00 1879 1486
## m1[1] 1.19 1.19 0.11 0.10 1.03 1.37 1.00 835 840
## m1[2] 4.33 4.32 0.29 0.27 3.88 4.82 1.00 881 874
## m1[3] 6.40 6.40 0.31 0.30 5.89 6.92 1.00 989 1039
## m1[4] 7.77 7.78 0.28 0.27 7.31 8.23 1.00 1262 1005
## m1[5] 8.68 8.69 0.25 0.24 8.27 9.08 1.00 1908 1583
## m1[6] 9.29 9.30 0.24 0.24 8.89 9.69 1.00 2367 1495
## m1[7] 9.70 9.70 0.26 0.26 9.27 10.12 1.00 1959 1465
## m1[8] 9.97 9.98 0.30 0.29 9.49 10.44 1.00 1483 1429
## m1[9] 10.15 10.15 0.33 0.33 9.61 10.69 1.00 1230 1259
## m1[10] 10.27 10.27 0.36 0.36 9.69 10.85 1.00 1107 1174
## y1[1] 1.19 1.17 0.91 0.87 -0.30 2.70 1.00 1758 2090
## y1[2] 4.32 4.33 0.95 0.93 2.78 5.85 1.00 1813 1812
## y1[3] 6.41 6.41 0.95 0.89 4.84 7.89 1.00 1926 2013
## y1[4] 7.76 7.76 0.94 0.91 6.24 9.33 1.00 1934 1947
## y1[5] 8.69 8.68 0.95 0.91 7.13 10.23 1.00 1899 1957
## y1[6] 9.27 9.25 0.96 0.91 7.76 10.91 1.00 1991 1825
## y1[7] 9.68 9.66 0.93 0.91 8.20 11.19 1.00 1945 1802
## y1[8] 9.98 10.01 0.98 0.91 8.38 11.60 1.00 2004 1916
## y1[9] 10.13 10.13 0.97 0.95 8.56 11.68 1.00 1935 1918
## y1[10] 10.28 10.27 1.04 1.03 8.61 11.99 1.00 1730 1932
ex8-5
sigmoid curve
y=Ym/ 1+exp(-b1* (x-b0)) -> y~B(Ym, 1+exp(-b1*(x-b0)))
b0,b1>0
x[0,Infinity), y[0,Ym]
(x=b0, y=Ym/2)data{
int N;
vector[N] x;
int Ym;
array[N] int y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0,upper=100> b0;
real<lower=0,upper=100> b1;
}
model{
y~binomial_logit(Ym,b1*(x-b0));
}
generated quantities{
array[N] int y0;
for(i in 1:N){
y0[i]=binomial_rng(Ym,inv_logit(b1*x[i]-b0*b1));
}
array[N1] int y1;
for(i in 1:N1){
y1[i]=binomial_rng(Ym,inv_logit(b1*x1[i]-b0*b1));
}
}n=20
x=runif(n,0,9)
ym=10
y=rbinom(n,ym,1/(1+exp(-2*(x-4))))
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,Ym=ym,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-5.stan')
mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.mcmc$metadata()$stan_variables## [1] "lp__" "b0" "b1" "y0" "y1"seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -37.93 -37.61 0.99 0.69 -39.92 -37.00 1.00 591 734
## b0 3.99 3.99 0.16 0.16 3.73 4.25 1.00 2032 1721
## b1 1.86 1.83 0.31 0.29 1.40 2.40 1.01 390 357
## y0[1] 0.05 0.00 0.22 0.00 0.00 0.00 1.00 1523 1503
## y0[2] 9.97 10.00 0.17 0.00 10.00 10.00 1.00 1629 NA
## y0[3] 8.17 8.00 1.32 1.48 6.00 10.00 1.00 1964 NA
## y0[4] 6.64 7.00 1.59 1.48 4.00 9.00 1.00 1981 1811
## y0[5] 9.24 9.00 0.87 1.48 8.00 10.00 1.00 1396 NA
## y0[6] 9.98 10.00 0.16 0.00 10.00 10.00 1.00 1998 NA
## y0[7] 0.28 0.00 0.54 0.00 0.00 1.00 1.00 1792 1756
## y0[8] 5.63 6.00 1.69 1.48 3.00 8.00 1.00 1862 1766
## y0[9] 0.13 0.00 0.37 0.00 0.00 1.00 1.00 1839 1898
## y0[10] 0.10 0.00 0.33 0.00 0.00 1.00 1.00 2070 2037
## y0[11] 9.97 10.00 0.17 0.00 10.00 10.00 1.00 2122 NA
## y0[12] 9.73 10.00 0.53 0.00 9.00 10.00 1.00 1976 NA
## y0[13] 0.47 0.00 0.72 0.00 0.00 2.00 1.00 1256 1383
## y0[14] 4.41 4.00 1.69 1.48 2.00 7.00 1.00 2101 2033
## y0[15] 9.97 10.00 0.16 0.00 10.00 10.00 1.00 1686 NA
## y0[16] 0.20 0.00 0.47 0.00 0.00 1.00 1.00 1760 1788
## y0[17] 9.99 10.00 0.12 0.00 10.00 10.00 1.00 2073 NA
## y0[18] 9.95 10.00 0.23 0.00 9.00 10.00 1.00 1896 NA
## y0[19] 8.97 9.00 1.04 1.48 7.00 10.00 1.00 1715 NA
## y0[20] 0.12 0.00 0.36 0.00 0.00 1.00 1.00 1720 1690
## y1[1] 0.05 0.00 0.23 0.00 0.00 0.00 1.00 1552 1560
## y1[2] 0.17 0.00 0.44 0.00 0.00 1.00 1.00 1546 1609
## y1[3] 0.61 0.00 0.83 0.00 0.00 2.00 1.00 1451 1528
## y1[4] 1.99 2.00 1.38 1.48 0.00 5.00 1.00 1650 1899
## y1[5] 5.07 5.00 1.72 1.48 2.00 8.00 1.00 1674 1753
## y1[6] 8.08 8.00 1.36 1.48 6.00 10.00 1.00 1893 NA
## y1[7] 9.42 10.00 0.77 0.00 8.00 10.00 1.00 1534 NA
## y1[8] 9.84 10.00 0.41 0.00 9.00 10.00 1.00 1696 NA
## y1[9] 9.95 10.00 0.22 0.00 10.00 10.00 1.00 1862 NA
## y1[10] 9.99 10.00 0.12 0.00 10.00 10.00 1.00 1795 NA
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,ymed=smy1$median,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=ymed),xy,col='black')
ex8-6
up and down
y=b0(exp(-b1* x)-exp(-b2* x)) -> y~N(b0(exp(-b1* x)-exp(-b2* x)),s)
b0,b1,b2>0, b1<b2
x[0,Infinity), 0<y<b0
(x=log b1-log b2 / b1-b2, y=max(y))data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0,upper=200> b0;
real<lower=0,upper=1> b1;
real<lower=0,upper=1> b2;
real<lower=0,upper=10> s;
}
model{
y~normal(b0*(exp(-b1*x)-exp(-b2*x)),s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0*(exp(-b1*x[i])-exp(-b2*x[i]));
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0*(exp(-b1*x1[i])-exp(-b2*x1[i]));
y1[i]=normal_rng(m1[i],s);
}
}n=20
x=runif(n,0,50)
y=rnorm(n,100*(exp(-0.03*x)-exp(-0.2*x)),1)
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-6.stan')
fn(mdl,data)## Initial log joint probability = -454.723
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 86 -6.55155 4.3167e-06 0.00416671 0.7339 0.7339 130
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -6.55
## b0 104.18
## b1 0.03
## b2 0.19
## s 0.84
## m0[1] 52.85
## m0[2] 30.98
## m0[3] 57.29
## m0[4] 40.25
## m0[5] 57.13
## m0[6] 36.59
## m0[7] 26.90
## m0[8] 51.95
## m0[9] 44.02
## m0[10] 25.50
## m0[11] 36.81
## m0[12] 52.08
## m0[13] 39.74
## m0[14] 35.14
## m0[15] 34.68
## m0[16] 55.62
## m0[17] 57.06
## m0[18] 58.09
## m0[19] 59.44
## m0[20] 44.26
## y0[1] 54.53
## y0[2] 31.69
## y0[3] 56.37
## y0[4] 40.12
## y0[5] 56.28
## y0[6] 37.88
## y0[7] 27.08
## y0[8] 50.76
## y0[9] 44.39
## y0[10] 24.05
## y0[11] 35.85
## y0[12] 53.21
## y0[13] 39.56
## y0[14] 34.84
## y0[15] 34.50
## y0[16] 55.36
## y0[17] 57.63
## y0[18] 57.86
## y0[19] 59.12
## y0[20] 44.67
## m1[1] 52.08
## m1[2] 60.07
## m1[3] 58.99
## m1[4] 54.45
## m1[5] 48.86
## m1[6] 43.26
## m1[7] 38.04
## m1[8] 33.34
## m1[9] 29.17
## m1[10] 25.50
## y1[1] 51.66
## y1[2] 60.22
## y1[3] 59.56
## y1[4] 53.04
## y1[5] 48.30
## y1[6] 42.19
## y1[7] 38.28
## y1[8] 32.44
## y1[9] 28.92
## y1[10] 23.97
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -10.87 -10.51 1.77 1.57 -14.43 -8.74 1.04 84 311
## b0 105.12 105.04 4.60 4.65 97.71 113.16 1.10 27 120
## b1 0.03 0.03 0.00 0.00 0.03 0.03 1.11 26 151
## b2 0.18 0.18 0.01 0.01 0.17 0.20 1.09 29 52
## s 1.04 1.02 0.20 0.17 0.76 1.43 1.14 23 82
## m0[1] 52.89 52.88 0.44 0.43 52.18 53.61 1.05 65 658
## m0[2] 30.93 30.92 0.40 0.39 30.27 31.60 1.06 52 497
## m0[3] 57.28 57.26 0.51 0.47 56.45 58.14 1.01 517 342
## m0[4] 40.24 40.24 0.29 0.27 39.76 40.72 1.00 1200 1034
## m0[5] 57.12 57.10 0.51 0.48 56.28 57.99 1.01 495 318
## m0[6] 36.56 36.56 0.31 0.29 36.03 37.06 1.02 944 713
## m0[7] 26.83 26.82 0.48 0.47 26.08 27.64 1.08 36 617
## m0[8] 51.99 51.98 0.43 0.42 51.31 52.68 1.04 69 570
## m0[9] 44.03 44.03 0.32 0.30 43.53 44.56 1.01 658 744
## m0[10] 25.42 25.42 0.50 0.49 24.63 26.29 1.08 34 658
## m0[11] 36.79 36.79 0.31 0.29 36.26 37.28 1.02 999 724
## m0[12] 52.06 52.04 0.63 0.58 51.03 53.17 1.03 120 196
## m0[13] 39.73 39.73 0.29 0.27 39.25 40.22 1.00 1232 901
## m0[14] 35.11 35.10 0.33 0.31 34.56 35.63 1.03 188 674
## m0[15] 34.65 34.65 0.34 0.31 34.10 35.18 1.03 144 674
## m0[16] 55.66 55.66 0.46 0.46 54.93 56.42 1.05 62 589
## m0[17] 57.06 57.04 0.51 0.48 56.21 57.93 1.01 488 318
## m0[18] 58.13 58.13 0.46 0.45 57.40 58.86 1.05 71 914
## m0[19] 59.48 59.48 0.44 0.43 58.77 60.20 1.04 96 1003
## m0[20] 44.27 44.27 0.32 0.31 43.76 44.80 1.01 626 730
## y0[1] 52.88 52.86 1.14 1.11 50.99 54.74 1.01 1453 1062
## y0[2] 30.94 30.93 1.13 1.10 29.08 32.77 1.00 1948 1884
## y0[3] 57.29 57.29 1.19 1.16 55.33 59.23 1.01 1151 1367
## y0[4] 40.25 40.21 1.09 1.02 38.39 42.05 1.00 1819 1556
## y0[5] 57.12 57.13 1.21 1.16 55.11 58.97 1.00 1616 1731
## y0[6] 36.53 36.53 1.10 1.10 34.68 38.29 1.00 1635 1466
## y0[7] 26.82 26.81 1.16 1.13 24.91 28.73 1.02 359 1590
## y0[8] 52.01 52.00 1.12 1.09 50.15 53.85 1.01 1183 1734
## y0[9] 44.05 44.03 1.08 1.03 42.26 45.86 1.00 2022 1929
## y0[10] 25.38 25.39 1.17 1.11 23.52 27.29 1.02 1109 1552
## y0[11] 36.78 36.81 1.11 1.09 34.92 38.54 1.01 2022 1752
## y0[12] 52.08 52.09 1.23 1.17 50.04 54.07 1.01 836 769
## y0[13] 39.76 39.78 1.11 1.07 37.91 41.59 1.00 2077 1810
## y0[14] 35.12 35.12 1.13 1.08 33.24 36.96 1.00 1878 1773
## y0[15] 34.66 34.65 1.15 1.11 32.85 36.58 1.00 1958 1664
## y0[16] 55.69 55.68 1.15 1.09 53.78 57.58 1.00 1407 1515
## y0[17] 57.02 57.04 1.16 1.11 55.19 58.94 1.00 1593 1650
## y0[18] 58.15 58.14 1.17 1.14 56.20 60.03 1.01 1417 1557
## y0[19] 59.46 59.45 1.13 1.08 57.63 61.33 1.00 1761 1804
## y0[20] 44.25 44.26 1.10 1.07 42.40 46.01 1.01 1296 1680
## m1[1] 52.06 52.04 0.63 0.58 51.03 53.17 1.03 120 196
## m1[2] 60.09 60.07 0.42 0.38 59.43 60.79 1.01 1170 1170
## m1[3] 59.03 59.04 0.45 0.44 58.32 59.76 1.04 83 1009
## m1[4] 54.49 54.48 0.45 0.45 53.76 55.23 1.05 62 518
## m1[5] 48.89 48.89 0.38 0.37 48.29 49.52 1.03 104 511
## m1[6] 43.26 43.26 0.31 0.29 42.77 43.78 1.00 758 828
## m1[7] 38.02 38.03 0.30 0.29 37.51 38.52 1.01 1191 745
## m1[8] 33.30 33.29 0.36 0.34 32.72 33.87 1.04 85 636
## m1[9] 29.11 29.10 0.43 0.42 28.41 29.85 1.07 42 499
## m1[10] 25.42 25.42 0.50 0.49 24.63 26.29 1.08 34 658
## y1[1] 52.02 52.05 1.20 1.15 50.04 53.98 1.01 1214 1736
## y1[2] 60.10 60.11 1.17 1.11 58.14 61.98 1.00 1774 1909
## y1[3] 59.03 59.03 1.12 1.12 57.21 60.84 1.01 1417 1640
## y1[4] 54.52 54.52 1.19 1.13 52.66 56.42 1.01 1313 1727
## y1[5] 48.88 48.86 1.14 1.07 47.00 50.65 1.00 1703 1932
## y1[6] 43.28 43.28 1.11 1.05 41.48 45.06 1.00 1886 1804
## y1[7] 38.02 38.02 1.08 1.05 36.17 39.80 1.00 1948 1690
## y1[8] 33.25 33.25 1.11 1.07 31.45 35.05 1.01 1678 1777
## y1[9] 29.12 29.14 1.12 1.08 27.28 31.01 1.01 1496 1733
## y1[10] 25.41 25.45 1.18 1.13 23.45 27.33 1.02 687 1133
